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TECHNICAL PUBLICATION 


STATISTICAL PROPERTIES OF MAXIMUM LIKELIHOOD ESTIMATORS 
OF POWER LAW SPECTRA INFORMATION 

1. INTRODUCTION 


A brief summary of the maximum likelihood estimation (MLE ) procedure developed for estimating 
power law spectral parameters in earlier works *•- begins with the probability density function of the 
astrophysics data set consisting of N detector responses v ( , e.g., energy deposit, as 


g(y,;0) = Jg(y ( - \E;p)<p(E;Q)dE , / = 1,--,jV, 

R 

where 0 denotes the vector of spectral parameters of an assumed energy spectrum <KE;Q) to be estimated; 
N is the number of detected events from observing range, R, of the instrument having response function, g, 
and energy resolution, p. Then the corresponding likelihood function is 


N 

L(0)=n 


1=1 


J g( \>j | £;p) <p(E\Q)dE 

R 


( 2 ) 


and the ML estimate of 0, say 0 ML , is chosen so that for any admissible value of 0, L(0 ml )> L(0) or 
equivalently, log[L(0 ML )] > log[L(0)]. 3 In practice, 0 ML can be obtained from equation (3) using an opti- 
mization routine, such as the Nelder-Meade simplex search algorithm, 4 to yield 


0 ML = min 0(0) , (3) 

{ 0 } 

where the objective function, O, in equation (3) is defined as 0(0) = -log[L(0)] for this minimization 
algorithm so that minimizing 0(0) maximizes log[L(0)] as desired, and where the integral appearing in 
equation (2) and thus intrinsic to equation (3) can be very accurately evaluated by numerical methods such 
as the method of Gaussian quadratures. 5 The Nelder-Mead algorithm does not require gradient information 
which is a vital consideration in selecting an optimization algorithm for this application because some 
energy spectra, such as the broken power law, are not differentiable everywhere. 


The MLE theory generally leads to lower bounds on the statistical errors (standard deviations) of 
the spectra information and the existence of such a bound, called the Cramer-Rao bound (CRB), is the 
bound below which the variance of an unbiased estimator cannot fall. This implies that irrespective of the 
method used to quantify the parameters from the data, there is a lower bound on the precision that cannot 
be superseded. 6 In the multiparameter case, if 0 is any unbiased estimator of 0, then 


var(0) > 


N 


9log[g(y;0)] 5log[g(y;0)]\l 1 

ae a© 7 ' /. 


(4) 


in the sense that the difference of these two matrices is positive semidefinite. The right-hand-side matrix in 
equation (4) is the CRB 7 and notationally will be referred to as I _l (0), where 1(0) is frequently called the 
information matrix; the notation < . > denotes “expected value”; and the superscript T stands for vector 
transpose. Thus, the variance of one component of is ft say ty, is bounded below as var(0 ; ) > 1^(0), where 
I” 1 (ft) is the /th diagonal element of I 1 (0) . When 0 consists of a single spectral parameter, e.g., a simple 
power law energy spectrum is assumed, the CRB is the right-hand side of the inequality 


var(0) > 


J( diog[g(y;fl)] y\ 

\l dO J / 


(5) 


Additionally, ML estimators generally possess the favorable large sample properties of consistency 
(unbiased) (PI ) and normality (P3). 6 - 8 - 9 This Technical Publication (TP) investigates the conditions whereby 
these two properties, along with efficiency which is attainment of the CRB and is referred to as property P2 
in this TP, are attained for an assumed simple power law energy distribution and a broken power law 
distribution, with emphasis on practical applications to instrument design and data analysis. 

1.1 Simple Power Law Energy Spectrum 

The simple power law suggests that the number of protons detected above an energy, E, is given by: 


N $ ( > E ) — M 



—(X j + 1 


( 6 ) 


where E is in units TeV, cc | is believed to be =2.8, and M A and E A are numbers associated with the detector’s 
collecting power (combination of size and observing time). In statistical terms, N s is assumed to represent 
an average number of events, while the actual number to be observed on any given mission would follow 
the Poisson probability distribution with mean number N s . The probability density function for galactic 
cosmic ray (GCR) event energy, E, is then given by 




<t> s (E) = 


(7) 


Ct] - 1 


,1-a. 


- E 


l-a. 


E~ a i 


for E\< E < Ei 


over an energy range [£j, £-,] that does not depend on the parameter a,. Because the actual incident 
particle energies are never observed, but only a measure of their energy deposition from their passage 
through the detector, the random variable Y is introduced to represent the detector's response, e.g., energy 
deposition, of a GCR proton of incident energy, E , and its stochastic response function, g, with energy 
resolution, p, which may or may not be energy independent. For specificity, a response function, g, based 
on simulation studies of a thin sampling calorimeter (TSC) concept for the Advanced Cosmic-Ray Compo- 
sition Experiment for the Space Station (ACCESS) will be used. It has a planned 3-yr program life cycle 
and is composed of a carbon target and sampling calorimeter. The TSC area is 1 m 2 with a target =0.7 
proton interaction lengths thick, sampled by X/Y pairs of square scintillating fibers. The fibers in the target 
are 2 mm thick and provide the approximate position of the interaction. The calorimeter consists of upper 
and lower parts totaling 25 radiation length (rl)-thick lead and contains 28 X/Y pairs of 500-p.m square 
scintillating fibers. The upper 3-rl-thick calorimeter is sampled each 0.5 rl, and the lower part is sampled 
each 1.0 rl. The total weight of the target and calorimeter is =2,600 kg, and the collecting power param- 
eters, M a and E a are estimated to be 160 and 500 TeV, respectively, implying that this TSC is expected to 
observe 160 proton events above 500 TeV over its expected 3-yr life cycle. 


The TSC performance predictions are based on the geometry and tracking particle physics simula- 
tion program (GEANT) simulations of energy deposition for monoenergetic protons at specified energies 
at 0. 1 , 1 , 10, 1 00, 1 ,000, and 5,000 TeV for this candidate detector. The Gaussian distribution was found to 
provide a reasonable description of the distribution of energy depositions at each of these incident ener- 
gies. 10 The mean detector response and the rms response were both found to be well approximated by a 
linear function of incident energy, E , in the range of interest for this study, which is typically in the range of 
20 to 5,500 TeV. Thus, the mean energy deposition, Y for a given incident energy, E, is defined to be 
= (a + bE) and the rms response defined as o^ E = ( c + dE), where the coefficients a, b, c, and d were 
estimated from the GEANT simulation results. 


Before investigating the properties of MLE for the TSC instrument, it is instructive to consider the 
concept of a zero-resolution instrument or so-called ideal detector because it sets an upper bound on the 
expected performance of any real detector of equal collecting power. This measurement bound is deter- 
mined by the CRB for the ideal detector which in turn establishes the limit in attainable precision with 
which unbiased spectra information can be obtained from a given science mission by any conceivable 
instrument with equal or less collecting power. Hence, it is useful in crafting realistic measurement goals 
for new science missions. Thus, the true character of the energy spectrum can be studied using an ideal 
detector without the additional statistical variations induced by real detectors, which is valuable when 
considering real detectors. 

An ideal detector's energy resolution, p, is equal to zero, so the standard deviation <5y\ E is zero for 
all GCR event energies, E. Hence, the incident GCR energies are precisely known from the inverse mean 
response, so that for the TSC having linear mean response gives E = (Y - a)/b , and using equation (5) 
provides the CRB as the right-hand side of the inequality: 


3 


var(d) > 


1 


N 


1 

^1+a, £ 1 +a, ( log[/ r 2] _ log[£j j) 2 T 

(Of, - 1) 2 

(E 2 E* x -E { E“') 2 [ 


( 8 ) 


and is asymptotically attained by the ML estimator. A key question then arises, “For what values of N is this 
asymptotic property P2, as well as PI and P3, achieved by MLE?” A battery of simulations was conducted 
to study this question consisting of 10,000 simulated missions for each of several values of A ranging from 
50 to 52,000* events per mission and with GCR energies from the energy range of 20 to 5,500 TeV. The 
ML estimate a ML was obtained for each mission by solving 


aiog[L] _ 1 

9cq cq - 1 


log[E 2 ]£ 


2 a| -log[E|]E| g| 
’ £.1— 


1 


N 


/=! 


(9) 


in terms of cq for the ideal detector and then also for a 40-percent resolution Gaussian detector (the TSC) 
by application of equation (3). For comparison, the estimation technique referred to as the “method of 
moments” is included for the ideal detector ( p = 0) for the case N = 200 events and N = 52,000 events, and 
consists of equating the sample mean E to the population mean as 


E = 


- 1 't 

6 ,- 2 ) 



( 10 ) 


and then solving the nonlinear equation (10) in terms of cq. 2 Figure 1 shows that MLE provides an unbi- 
ased estimate of cq when N > 1 ,000, but with an ever-increasing bias as the number of events diminishes. 
Note that even though a ML is biased when N = 200 for the 40-percent resolution Gaussian TSC, its bias is 
significantly less than the bias of the method of moments estimator cq for the ideal detector having perfect 
energy resolution. 


An analytical expression that allows one to compute the bias of a ML for the ideal detector can be 

constructed by noting in equation (9) that « ML is a function of the logarithm of the event energies by the 

1 ‘ v 

term TrX'osfc] . Thus, the random variable IT = log[£] is introduced having probability density 

TV /=| 

function 


/(»’> 


(cq -l)£, ai £“V v(1 ~ ai ) 

fcj | lL 'y L j h. ^ 


log[E ( ]< w< log[E 2 ] 


( 11 ) 


*The TSC used in the ACCESS concept study would detect 52,000 events on average over the energy 
range of 20 to 5,500 TeV when the spectral parameter, cq, is assumed to be 2.8. 
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♦ Ideal Detector □ 40% Resolution, Gaussian 

▲ Method of Moments (Ideal) Theoretical Mean (Ideal Detector) 



Figure 1 . Mean estimates of a ML as a function of the number of events, N, 

showing a bias for N<\ ,000. Triangle marker at /V = 200 and 52,000 
is for method of moments. 


and with mean and variance 




= log[£|J + 


1 

— 1 + oc j 





<5~ w 


1 

(~l + oc l)~ 



E^eU 2 

or 


( 12 ) 


1 N 

where S = log[Ej]- logffT], co = -E® 1 E-, + £, £?' , and, by the central limit theorem, — X w / is nor ' 

1 ; ' =1 

mally distributed with mean jt%and variance — c 7^ . Consequently, the probability distribution of the ML 
estimator a ML of a { using the ideal detector is obtained by solving the probability equation. 


1 

log[£ 2 ] £ 2 -aML -logt^]]^ 1 aMl 

-Vw 

' 

£|“ a ML _ 

a ML ~ 1 


°w 



4n 




— co 


( 13 ) 
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in terms of a ML for various values of Z. Letting Z vary from -5 to 5, setting N to each of the number of 
events N used in figure 1, and then numerically evaluating the mean of a ML from the probability density 
function constructed from equation (13) gives the solid curve shown in figure 1 , indicating good agreement 
with the simulation results. 

An interesting observation from figure 1 is that one would likely conclude, and incorrectly, is that a 
significant difference between the slopes of two cosmic-ray elemental species exists if their respective 
number of events were significantly different from each other and at least one had fewer than 1 ,000 events. 
This is because for two given cosmic-ray elemental species, A and B, with simple power law parameters, a 
and ft, the hypothesis H 0 : a- (3 = 0 (same “slopes”) versus H { :a- /?*0 uses the test statistic 
that will inherit the bias(s) shown in figure 1 when N is < 1 ,000. Thus, an interesting study would be to plot 
the estimate of the slope parameter for each of several elemental species as a function of A comprising their 
respective data sets to see if it resembles figure 1 . It should also be understood that this bias as a function of 
N would be even worse had the method of moments estimation procedure been used to estimate the spectral 
parameters as previously noted in figure 1 . 

A comparison of the standard deviation of a ML for the ideal detector with the CRB determined 
from equation (8) as a function of the number of events A is depicted in figure 2 and clearly shows that a ML 
attains the CRB for A > 100. The standard deviation of the method of moments estimator cq for the ideal 
detector is also provided for comparison for A = 200 and A = 52,000, and note once again that MLE 
provides a superior estimator. 


Standard Deviation ot ML Estimator and CRB for Ideal Detector 

— Standard Deviation (Ideal) CRB a Method of Moments (Ideal) 



Number of Events 


Figure 2. Standard deviation of a ML and the CRB as a function of A. 

The comparison for A = 52,000 (rightmost markers in fig. 2) is somewhat visually misleading 
because of the scale of the vertical axis. Their actual values are 0.008, 0.012, and 0.008 for a ML , cq , and 
the CRB, respectively, and the ratio 0.012 to 0.008 is 1.5 for the method of moments and 1.0 for MLE, 
which is a measure of the efficiency of the estimators and again shows MLE is far superior to the method 
of moments. Simply put, the improvement in measurement precision provided by MLE over the method of 
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moments can be roughly equated to doubling the collecting power of the instrument, because doubling the 
collecting power reduces the standard deviation by 1 / y2 when the CRB is attained. Furthermore, 
this ratio of ~1.5 remained steady as the detector resolution varied from zero to 50 percent, as shown in 
figure 3. This fact, coupled with the better performance in achieving PI as previously discussed, explains 
why MLE is superior to the method of moments when estimating power law spectra information. 

ML Estimator a ML , Method ot Moments Estimator, 
and CRB Versus Detector Energy Resolution. 

M=52,000 Events 

♦ Method of Moments ■ Maximum Likelihood CRB 



Figure 3. MLE and method of moments as a function of instrument energy resolution 
for the proposed TSC with Gaussian response function. 


Property P3 is investigated using a frequency histogram of a ML based on the 10,000 simulated 
missions when N = 50 events per mission and shows a significant right-hand skewness (fig. 4a), and thus, 
a clear departure from normality (Gaussian fit is illustrated as smooth curves in fig. 4), while a similar 
comparison for the case N = 52,000 shows is very normally distributed. Visual studies of the interme- 
diate values of N showed the frequency histograms to be normally distributed in appearance for N> 1 ,000 and 
is in concert with the bias study depicted in figure 1 . 


Histogram of ML Estimate of a, With Gaussian Fit, 
Gaussian Detector, N= 50 



2.5 3 3.5 4 

ML Estimate of a, 


4.5 


Histogram of ML Estimate of a t With Gaussian Fit, 
Gaussian Detector, N= 52,000 


(at 



Figure 4. Frequency histograms of f° r 10,000 simulated missions with (a) N - 50 

and (b) N = 52,000 events, using the TSC with its 40-percent resolution Gaussian 
response function. 
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A relative frequency histogram of a ML based on 10,000 simulated missions with N = 50 per mis- 
sion using an ideal detector having zero energy resolution is shown in figure 5. Also shown is the theoreti- 
cal distribution of a ML obtained from equation (13) with parameters set to N = 50, £j = 20 TeV, £\ = 5,500 
TeV, and a j = 2.8 and illustrates the close agreement between simulation and theory. 

Histogram of ML Estimate of a , With Theoretical Distribution When N = 50 for the Ideal Detector 



Figure 5. Frequency histogram of or ML based on 10,000 simulated missions with N = 50 
using an ideal detector. The smooth curve is the theoretical distribution of a M1 
obtained from equation (13). 

Solving equation (3) to obtain the ML estimate for the case where the events are measured by a real 
detector having nonzero energy resolution is straightforward, and checking consistency (PI ) and normality 
(P3) is easily performed. However, checking efficiency (P2) can be quite formidable because of the term 


/ ( dlog[#(y;«i )]' 

~\ 7 

a , 

’e 2 

J g(y \E,p)<f>s(E;a i )dE 

^•2 

r -A 

E, 

r 

\ dct\ 

H 

3o, 108 


J g(y | E,p)<p s (E;a l )dE 

\ V i ) 

/ 0 

V 

- £ ' 


Ei 

) 


required to compute the CRB, coupled with the fact that the detector response function, g, in equation ( 14) 
can be quite complicated. For example, g could be Gaussian with non-negativity constraint y > 0 and 
energy-dependent resolution function, p(E ), that in turn requires an energy-dependent normalizing coeffi- 
cient. Fortunately, equation (14) can be numerically evaluated using the symmetrized form of the numeri- 
cal derivative, 1 1 


r(l) _ /<■* + *>-/<*-*> ,15) 

J ' 2 h 

to approximate the derivative in equation (14) and in conjunction with the method of Gaussian quadratures 
to calculate the definite integrals. The fact that the CRB in equation (8) for the ideal detector must match 
the CRB obtained from equation ( 14) when the detector resolution is zero (p—> 0) provides a means to tune 
the numerical differentiation parameter, h, and the integration parameters, e.g., the upper integration limit 
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for y as well as the number of partitions used in the numerical integration in both integration variables, E 
and y, in equation ( 14). For the TSC instrument with a data analysis range of 20 to 5,500 TeV, setting h to 
0.0001 and the upper integration limit of y to 35,000 GeV in place of infinity in equation (14), and using 
10-point Gaussian quadratures over subintervals over both integration ranges provides the somewhat sur- 
prising result of 13-decimal-place accuracy in the numerical evaluation of equation (14) when compared to 
the exact value obtained from equation (8) for the ideal detector. This accuracy in the numerical evaluation 
of equation ( 14) was independently confirmed using the numerical integration routine in MATHEMATICA® . 

Figure 6 illustrates the convergence of the standard deviation of the ML estimator to the CRB 
computed using equation ( 14) for a 40-percent resolution Gaussian detector as a function of N. The stan- 
dard deviation of a ML is based on a battery of 10,000 simulated missions for each value of N, where N 
ranges from 50 to 52,000 events per mission. The CRB for the ideal detector is included as a reference 
curve ( — ). 


• Standard Deviation (40% Gaussian) — CRB (40% Gaussian) — CRB (Ideal) 



Figure 6. Standard deviation of the ML estimator versus N for a 40-percent resolution 
Gaussian detector based on 10,000 simulated missions at each value of N and the 

CRB for the 40-percent Gaussian detector computed using equation ( 14) ( ). 

The CRB for the ideal detector is included as a reference curve ( ). 

When MLE is being used in the design phase of an instrument to estimate its expected performance 
and if the simulations indicate that MLE does in fact provide unbiased spectra information and approxi- 
mate attainment of the CRB for the science mission under study, then equation ( 14) can be used to evaluate 
the relative merits of various instrument design parameters without performing additional simulations. 
This has tremendous practical value in design parameter trade studies because equation (14) can be evalu- 
ated in mere seconds, while the equivalent information from Monte Carlo simulations can take several 
days. For example, because we know (from Fig. 2) that MLE attains the CRB for N > 100 events, equation 
(14) can be used to compute the family of curves shown in figure 7 that relate the precision with which a, 
can be measured as a function of detector resolution and collecting power. This implies instrument design- 
ers should first attempt to maximize collecting power and then improve resolution, and in that order. The 
proposed TSC instrument, with its expected 52,000 events, is indicated by the square in figure 7. The reader is 
referred to figure 3 for a detailed view of the CRB calculated using equation ( 14) for the TSC as a function 
of detector energy resolution. 
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CRB for a 1 Using Gaussian Detector, 0%-40% Resolution 



Number of Events 


Figure 7. CRB as a function of N for detector energy resolutions in the range 0 < p <0.40. 

The proposed TSC instrument with its expected 52,000 events is indicated by the 
square, and the observing energy range was set to 20-5,500 TeV. 


A detailed simulation study of the TSC-sized ideal detector, with its expected N = 52,000 events for 
the observing range 20-5,500 TeV and with a, = 2.8, was conducted and a ML obtained from equation (9) 
for each of 1 million missions (each mission detected 52,000 events), yielding a mean and standard devia- 
tion value of a ML to be 2.80003 and 0.007905, respectively. Constructing the probability density function 
of a ML from equation (13) and then numerically evaluating its mean and standard deviation gives 2.80003 
and 0.00791 1, respectively, while the CRB calculated from equation (8) gives 0.00790998, illustrating the 
remarkable agreement between simulation and theory and attainment of P2. Note that a ML is essentially 
unbiased too and thus, PI is approximately attained. Because 5.2 x 10 10 random numbers were required for 
this 1 -million simulated missions study, it is crucial to use a random number generator having a period 
longer than 5.2 x 10 10 , such as the generator used in this study which has a period of = 10 18 . 

1.2 Broken Power Law Energy Spectrum 

The broken power law energy spectrum suggests a transition from spectral index a, below the knee 
location at energy E k to a steeper spectral index > a x above the knee.* The broken power law spectrum 
predicts that the number of protons detected above an energy, E, is given by: 


*The case a 2 <oc j and where E k is referred to as the ankle can also be handled by this MLE procedure. 
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N b (>E) = \ 
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\ e a 


V -a i + 1 


{ 17 \- a 2 + 1 


yr-kJ 


for E > E^ 


\ a 2~ l J 

N s (> E)-[N S (> E k )-N B (> E k )] for E < E k 


(16) 


where E is in units TeV, M A and E A are 160 and 500 TeV, respectively, as before for the TSC instrument, 
N S (>E) is defined in equation (6), and currently available measurements suggest that a, is =2.8, Cb is 
thought to be somewhere between 3. 1 and 3.3, and E k is parameterized in the range of 100 to 300 TeV for 
this study. The broken power law probability density function (j) B is obtained by normalizing N B over an 
observing range [£j £3>] of interest and is defined in equation (21 ) of appendix A. 

The likelihood function of a random sample of N Galactic Cosmic Rays (GCR) events from the 
broken power law spectrum detected by the ideal detector having perfect energy resolution, regarded as a 
function of the unknown vector of parameters 0 = (0f|, o^, E k ), is 


L(0) = A(0) 


N 


nf 


\ E i <E k 


-a i 


n ^ 


\-cc 2 


K E j* E k 


Ei, 


E { < E h Ej < E 2 


(17) 


where the first product is over the event energies below the knee location E k and the second product is over 
those event energies above E k , and they total in number to N , and A(0) is the normalizing coefficient given 
in equation (22). 

The Nelder-Mead simplex method can then be used to obtain 0 ML from equation (18), where the 
objective function 0(0) is defined as minus the log-likelihood function, so that 


©ml = mini 

e 


-iVlogA(0) + cq 


X lo § 


E i< E k 


It 

E l . 


+ a-> 


X '°g 


E J 


\ E j- E k 


(18) 


for N events detected by the ideal detector, while equation (2) must be used to construct the likelihood 
function for a real detector having response function, g, and energy resolution, p. with 1 V instrument 
responses y, . Consequently, the ML estimate ®ML * s 


(V 

/ 

' E k 

#ML = '"“'X 

-log 

f g(yi 1 E >p) M0) 

(•)5 


Li. 1 


where the range of integration must be split at E k 
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J g{}’i I £,p) A(0) — dE 
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(19) 


at each step in the simplex search for 0 ML = (a,, a 2 , E k ) ML - 


1 1 



To numerically explore the properties of 0 ML for the broken power law distribution, the vector of 
spectral parameters is first set to 0 = (2.8, 3.3. 100 TeV) and events simulated from the energy range 20 to 
5,500 TeV for each of several values of N selected so as to provide an average of 500, 1,000, ..., 5,000 
events above E k as shown in table 1 . The notation N 2 is introduced to denote N B {>E k ), N for N B (>E l ) so 
that N\=N-N 2 , and the notation N(N 2 ) means “a total of N events, of which N-, of them are above the 
spectral knee E k ." 

Table I . Number of events used in broken power law simulations for 0.5 break-size study. 



10,939 

21,877 

32,816 

43,754 

65,631 

87,508 

109,390 

n 2 

500 

1,000 

1,500 

2,000 

3,000 

4,000 

5,000 

N 

11,439 

22,877 

34,316 

45,754 

68,631 

91,508 

114,390 


For each value of N in table 1, 10,000 missions were simulated and for each of these missions, 0 ML 
was obtained using equation (18) for an ideal detector and equation (19) for the TSC detector having 
Gaussian response function g and constant 40-percent energy resolution over the simulated energy range of 
20 to 5,500 TeV. 

Figure 8 depicts the mean of the ML estimates of a { and Ob versus the number of events N used in 
the simulations and shows that when the collecting power of the detector provides >1,500 events above E k 
(corresponding to third set of markers from left), property PI is essentially attained by the TSC instrument 
since the relative bias is <3 percent for the 40-percent resolution Gaussian detector, and is even better for 
the ideal detector having zero energy resolution. 

Similarly, figure 9 illustrates the bias (recall E k = 100 TeV in these simuations ) of the ML estimate 
of E k as a function of N for the TSC instrument and the ideal detector. Note that property PI is again 
roughly attained by the TSC (relative bias is <2 percent) when there are >1,500 events above 

♦ «, (Ideal) ■ a 2 (Ideal) 

a, (40% Gaussian) — #- a 2 ( 40 % Gaussian) 
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E k (Ideal) E k (40% Gaussian) 



Figure 9. ML estimate of E k as a function of detector collecting power using a 40-percent 
resolution Gaussian response function (the TSC) and the ideal detector. 

The knee location E k was set to 100 TeV in these simulations. 

Next, property P2 is investigated and requires the construction of the 3-by-3 information matrix 
1(0). Equation (32) of appendix A provides 1(0) for the ideal detector, while for a real detector with 
response function, g, and energy resolution, p, the //'-element of 1(0) is, by equation (4), 




d i 

m' og 
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jg(y\E,p)<t> B (E;Q)dE 

E\ 

| g(y | E,p) <p B (E;Q)dE 


\dv 


d 

*£> 

— log 

dd & 

j g{y\E,p)<p B (E\Q)dE 

- E > J. 


( 20 ) 


and can be accurately computed using the numerical methods discussed in the simple power law section, 
and where the notation in equation (20) defines 0, = a x , 0 2 = « 2 ' and e ?,= E b and where the integration range 
[£,, E-,] must be split at E k for the inner three integrals. 

A comparison of the CRB obtained from equation (20) for (X\ using the TSC with its 40-percent 
Gaussian response function with the simulation results is shown in figure 10. Note the CRB is attained 
when the number of events above the knee location is >1,500. The case 1 1,439 (500) had several simulated 
missions in which the MLE procedure gave an estimate 0 ML of 0 that suggests a simple power law would 
probably be an adequate explanation of the simulated events. These Errant Estimates were characterized as 
EE1 ) E k and a s are both very large relative to their assumed values of 100 TeV and 3.3 in the simulations, 
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and EE2; E k and a, are both very small. While the condition a } ~ is normally associated with suggest- 
ing a simple power law adequately fits the data, these unlucky missions illustrate the beauty of the MLE 
procedure in finding two other conditions whereby a broken power law collapses into a simple power law. 
The first condition is a broken power law with E k above the range of detected events and a-,— »°o in an effort 
to explain the absence of events above E k , which is indeed just a simple power law over the range of 
detected events and implied by EE1. The second condition is a broken power law with E k below the range 
of detected events and a } — >0 and implied by EE2. Eliminating these errant estimates of a, gives the 
trimmed standard deviation depicted at A = 1 1 ,439 (and A 2 = 500) in figure 10 and symbolized by a filled 
circle on the plot. The CRB for the ideal detector calculated from equation (42) is also shown in figure 10 
with corresponding simulation results. Additionally, the difference between the covariance matrix of the 
ML estimates and I~'(0) was noted to be positive definite as each its three eigenvalues were positive, with 
two of them approximately zero for all values of Ain table 1 used in the simulations. 

Similar results are illustrated in figure 1 1 for and figure 12 for E k . Trimmed estimates for these 
two make little difference because of their already larger variance relative to that of a j . 


Standard Deviation of ML Estimate of a, and CRB 


ao (Ideal) CRB (Ideal) •Trimmed 

■ o (40% Gaussian) — — CRB (40% Gaussian) 



Total Number of Events N 


Figure 10. CRB of a { using TSC ( ) and ideal detector ( (obtained from equation (20) 

versus collecting power. Standard deviations of ML estimates from simulations for 
values of A in table 1 indicated by markers. 


14 




Standard Deviation of ML Estimate of a 2 and CRB 

■ o (40% Gaussian) CRB (40% Gaussian) 

ao (I deal) CRB (Ideal) •Trimmed 



Figure 11. CRB of o ^ using TSC ( ) and ideal detector ( ) obtained from equation (20) 

versus collecting power. Standard deviations of ML estimates from simulations for 
values of N in table 1 indicated by markers. 


Standard Deviation of ML Estimate of f^and CRB 

a o k (40% Gaussian) CRB (40% Gaussian) a o k (Ideal) CRB (Ideal) 



Figure 12. CRB of E k using TSC ( ) and ideal detector ( ) obtained from equation (20) 

versus collecting power. Standard deviations of ML estimates from simulations for 
values of N in table 1 indicated by markers. 
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The property of asymptotic normality (P3) of 0 ML is next investigated with the aid of relative 
frequency histograms of the components of 0 ML provided from the simulations. Figure 1 3 shows relative 
frequency histograms of the 10,000 ML estimates of a, and a 2 for the two collecting powers that provide 
1 1 ,439 (500) events and also 1 14,390 (5,000) events and correspond to the first and last column of table 1 . 
As before, the detector here is the TSC with its Gaussian response function and 40-percent energy resolu- 
tion and with its collecting power adjusted through the choice of N. Note that while the histograms corre- 
sponding to the larger collecting power are approximately normally distributed and well separated, those 
corresponding to the smaller detector are skewed and even slightly overlapping, indicating the onset of 
difficulties in detecting the broken power law parameters. Relative frequency histograms for the ideal 
detector (not shown) show no overlap for the N= 1 1 ,439 (500) case and suggest that this is the approximate 
boundary for fixed N where detector resolution can play a leveraging role for this set of parameters. 


Histograms of ML Estimate of a, and ctj, for Gaussian Response Function With 40% Resolution. 
6 = (2.8, 3.3, 100 TeV), Energy Range 20-5,500 TeV. Average of 500 and 5,000 Events 
Above Knee, With 11,000 and 110,000 Below, Respectively. 

10,000 Simulated Missions 



Figure 13. Relative frequency histograms of the ML estimate of a x (leftmost two histograms) 
for N = 11 ,439 (broadest of the two) and N = 1 14,390 (narrow histogram). 

Rightmost two histograms similarly defined for a 

Frequency histograms of the ML estimates of E k for these two cases of N are shown in figure 1 4 and 
once again, note that the larger sized detector has roughly attained P3 while the smaller sized detector has 
not, and in fact a "bump" in the tail of the broader distribution for the smaller detector is seen, suggesting 
a simple power law would likely be an adequate explanation of these particular mission results. This 
situation was previously discussed and referred to as EEL 
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Histograms of ML Estimate of for Gaussian Response Function With 40% Resolution. 
0 = (2.8, 3.3, 100 TeV), Energy Range 20-5,500 TeV. Average of 500 and 5,000 Events 
Above Knee, With 11,000 and 110,000 Below, Respectively. 

10,000 Simulated Missions 



Figure 14. Relative frequency histograms of the ML estimates of E k for A/ — II ,439 
(broadest of the two and with bump in right-hand tail) and N = 1 14,390. 

Next, figure 1 5 shows relative frequency histograms of the 10,000 ML estimates of oq and 0^ when 
N = 22,877, providing an average of N 2 = 1,000 events above E k , and the two histograms are seen to be 
clearly separated. This suggests that a detector with this collecting power and a 40-percent resolution 
Gaussian response function could indeed measure the three broken power law spectral parameters when 
their true values are 0 = (2.8, 3.3, 1 00 TeV). Because the concept TSC that was studied would detect 5 1 .576 
(2,255) events on average in the energy range 20-5,500 TeV when 0 = (2.8, 3.3, 100 TeV), it is concluded 
that it could measure the three spectral parameters when E k ~ 100 TeV and the break-size is =0.5. 


Histogram of ML Estimate ot a, and a 2 for Gaussian Response Function With 40% Resolution. 
0 = (2.8, 3.3, 100 TeV), Energy Range 20-5,500 TeV. Average of 1,000 Events 
Above Knee, 21,877 Below. 

10,000 Simulated Missions 



2.5 2.8 3.1 3.4 3.7 4.0 


Spectral Parameters a, and oq 

Figure 15. Relative frequency histograms of ML estimates of Gq and G^> for N = 22.877 and N 2 = 1,000. 
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1.3 Break-Size 0.3 Study 


For this study the vector of spectral parameters is set to 0 = (2.8, 3.1, 100 TeV) and events 
simulated from the energy range 20-5,500 TeV for each of several values of N selected so as to provide an 
average of 1 ,000, 2,000, .... 5,000 events above E k as shown in table 2. (Values of N-, < 1 ,000 produced too 
many errant ML estimates 0 ML of 0 to be useful.) 

For each value of N in table 2, 5,000 missions were simulated and for each mission, 0 ML 
was obtained using equation (18) for an ideal detector, and equation (19) for the TSC detector having 
Gaussian response function and constant 40-percent energy resolution over the simulated energy range 
20-5,500 TeV. Figure 1 6 shows that when the number of detected events above the knee is >2,000, the ML 
estimate of a, and a 2 is essentially unbiased and property PI is attained, while figure 17 indicates the ML 
estimate of E k is still somewhat biased, even when N 2 = 2,000 (second markers form left) for the 
40-percent resolution Gaussian detector, which is perhaps not surprising in light of the more difficult 
estimation task for this smaller break-size case. 


Table 2. Number of events used in broken power law simulations for 0.3 break-size study. 



19,977 39,954 59,931 79,909 99,886 


1.000 2,000 3,000 4,000 5,000 

N 

20,977 41,954 62,931 83,909 104,886 


Average ML Estimate oi a, and a 2 

— ♦- cq (Ideal) — a- a z (Ideal) a, (40% Gaussian) a 2 (40% Gaussian) 



Figure 16. ML estimate of a, and as a function of detector collecting power 
when the spectral break size is 0.3 for the TSC and the ideal detector. 
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Average ML Estimate of E k 

—tr~ E fr (Ideal) E k (40% Gaussian) ♦ TSC (proposed size) 



Figure 17. ML estimate of E k as a function of detector collecting power using the TSC 
and ideal detector when the break size is 0.3. The concept TSC with its 
expected A = 5 1 ,790 events is indicated by the diamond and is based on 25,000 
simulated missions (others are 5,000 missions each) and suggests the marker to its 
immediate right is probably a little on the high side. 

Figure 1 8 shows the CRB of a, using the TSC detector and also the ideal detector versus the 
number of detected events N. The markers represent the standard deviation of the 5,000 ML estimates of 
a, based on the simulations. Note that when N 2 = 1,000, MLE experienced several missions resulting in 
errant estimates of a ( in its attempt to place the knee before the data and then drive a { — >0 (condition EE2), 
suggesting a simple power law might be a suitable fit for those simulated missions. Trimmed estimates are 
also provided in figure 18 corresponding to the cases where N 2 = 1,000 and N 2 = 2,000 and indicated by 
filled circles. Also note the ideal detector with its zero-percent resolution attains the CRB for all the values 
of N in table 2. 

Similar comparisons between the CRB and the ML estimate of Ch and E k are shown in figures 19 
and 20, respectively. Figure 20 indicates the CRB for E k is particularly difficult to attain, even for the ideal 
detector. Trimmed estimates are indicated by filled circles for the first two values of A, corresponding to an 
average of 1 ,000 and 2,000 events above E k , respectively. 
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Standard Deviation 


Standard Deviation of ML Estimate of a 1 and CRB 

o (Ideal) CRB (Ideal) • 0.5% Trimmed 

■ o (40% Gaussian) CRB (40% Gaussian) ♦ TSC (Proposed Size) 



Figure 1 8. CRB of a i using TSC (— — ) and ideal detector ( ) 

versus collecting power. Standard deviations of ML estimates 
from simulations for values of N in table 2 indicated by markers. 


Standard Deviation of ML Estimate of a 2 and CRB 

■ g (40% Gaussian) CRB (40% Gaussian) A o (Ideal) 

• 0.5% Trimmed CRB (Ideal) « TSC (Proposed Size) 



Total Number of Events N 


Figure 19. CRB of a 2 using TSC ( ) and ideal detector ( ) 

versus collecting power. Standard deviations of ML estimates 
from simulations for values of N in table 2 indicated by markers. 
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Standard Deviation of ML Estimate of E k and CRB 


a o k (40% Gaussian) CRB (40% Gaussian) * o*(ldeal) 

• 0.5% Trimmed CRB (Ideal) ♦ TSC (Proposed Size) 



Figure 20. CRB of E k using TSC ( ) and ideal detector ( ) 

versus collecting power. Standard deviations of ML estimates 
from simulations for values of N in table 2 indicated by markers. 

To investigate the properties of 0 ML for the proposed TSC with its 40-percent resolution Gaussian 
response function, 25,000 missions were simulated with 0 = (2.8, 3. 1 , 100 TeV), providing 5 1 ,790 (2,470) 
events on average from the observing range 20-5,500 TeV. Frequency histograms of the ML estimates of 
a, and Oh are shown in figure 21 for the proposed TSC and a clear separation between the histograms is 
observed. A slight right-hand skewness in the ML estimates of ts noted. 

Relative Frequency Histogram of ML Estimate of a, and ctj 
e = (2.8, 3.1 , 100 TeV), Observing Range: 20-5,500 TeV, 

Gaussian Response Function With 40% Resolution. 

25,000 Missions 



Figure 21 . Relative frequency histograms of ML estimates of a, and a 2 for the proposed TSC. 
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Fi S ure 22 shows the relative frequency histogram of the ML estimates of E k using the proposed 
TSC and the long, right-hand tail suggests the possibility of a few missions that resulted in errant estimates 
of the form EEL Also note the skewness and thus slight departure from normality (property P3). 

Table 3 gives summary statistics of 0^^ h)r the simulated missions. The rows labeled "theoretical 
limits” under the Comments column provide the input parameters 0 used in the simulations and the CRB 
obtained from equation (20), indicating that 0^ is approximately unbiased, efficient, and normally dis- 
tributed so that properties PI, P2, and P3 are roughly attained by the proposed TSC for this set of spectral 
parameters. Similar information for the zero-resolution ideal detector is also provided in table 3. 
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Figure 22. Relative frequency histograms of the ML estimates of E k for the proposed TSC. 


Relative Frequency Histogram of ML Estimate of E k 
9 = (2.8, 3.1 , 100 TeV), Observing Range: 20-5,500 TeV, 
Gaussian Response Function With 40% Resolution. 
25,000 Missions 



Table 3. Summary statistics based on 25,000 simulated missions with 0 = (2.8, 3.1, 100 TeV) 
and observing range 20-5,500 TeV, for the TSC having Gaussian response function. 


Resolution 

(%) 

f* (TeV) 

N(NJ 

Mean 

Standard Deviation 

Comments 


«2 

£*(TeV) 


a 2 

CfeV) 

0 

100 

51.790 (2,470) 

2.80 

3.10 

100 

0.012 

0.043 

10.6 

Theoretical limit 

0 


2.80 

3.10 

101 

0.012 

0.048 

13.7 

Simulation 
(25,000 missions) 

40 

100 

2.80 

3.10 

100 

0.020 

0.061 

21.2 

Theoretical limit 

40 

51,790(2,470) 

2.80 

3.11 

103 

0.022 

0.068 

24.3 

Simulation 
(25,000 missions) 





The difficulty of the MLE task is only partly appreciated from the preceding study of the attainment 
of (or lack of) the three statistical properties — PI, P2, and P3. Figure 23 illustrates the objective function 
in the vicinity of 0 ML for the ideal detector for a particular simulation mission in which 0 = (2.8, 3.3, 
100 TeV) and there were 1 14,385 events from the energy range 20-5,500 TeV of which 4,819 were above 
the knee at 100 TeV, with equation (18) yielding 0 ML = (2.805, 3.319, 95.16 TeV). Figure 23 shows the 
objective function in the neighborhood of 0 M l- keeping a x fixed at 2.805 and letting and Ej. vary in the 
region around 0 ML . Note the surface is well behaved and the minimum is easily found (and hence, 0 ML ), 
and is representative of all the surfaces that were viewed when properties P 1 , P2, and P3 are approximately 

attained. 


Objective Function Versus a 2 and f*in Vicinity ol e ML , With oq Fixed at its ML Estimate 2.805 



Figure 23. Objective function given by equation (18) for the ideal detector in the neighborhood 
of 0 ML for a simulated mission, keeping a x fixed at its ML estimate. There were 
1 14,385 (4,819) events in this mission. The vertical axis has been scaled to 1. 
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The alignment or tilt of the surface in figure 23 is interesting and the contour plot in figure 24 
illustrates the strong correlation between os and E k which is a direct consequence of the mathematical 
formulation of </> g in equation (2 1 ) where the knee E k acts as a “hinge” connecting the two distinct parts of 
the spectrum, and one can easily visualize a correlation between ctj and E k as well, whereas ctj and cc~, 
appear to be only slightly correlated (not shown). The information matrix in equation (32) for the ideal 
detector can be used to approximate the correlations for this case and gives zero between a, and a,. 0.41 
between otj and E k , and 0.64 between os and E k . 


Contour Plot of Objective Function Versus a 2 and E k in a Neighborhood 
of 0 ML . With a, Fixed at its ML Estimate 
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Figure 24. Contour plot of the objective function for ideal detector in the neighborhood 
°f ®ml for a simulated mission, keeping a ] fixed at its ML estimate. There 
were 1 14,385 (4,819) events in this mission. 
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On the other hand, considering a much smaller ideal detector that detected 4,575 events of which 
207 were above the knee on a particular simulated mission, equation (18) yields 0 ML = (2.856, 3.354, 
156.57 TeV) and corresponds to the broader of the two minima shown in figure 25 which reveals a much 
more difficult estimation terrain in the neighborhood of 0 ML and in fact, suggests the possibility of mul- 
tiple optimal solutions and confidence intervals for and Ej. that would necessarily consist of disjoint 
subintervals! 


Objective Function Versus a, and in Vicinity of 0 ML , With a, Fixed at its ML Estimate 



Figure 25. Objective function for ideal detector in the neighborhood of 0 M ^ f° r a simulated 

mission, keeping a, fixed at its ML estimate. This mission consisted of 4,575 (207) 
events. The vertical axis has been scaled to 1 . 
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As observed in figures 13 and 15, when histograms of the ML estimates of and begin to 
overlap, it clearly signals the onset of difficulties in estimating the broken power law spectral parameters 
using MLE (and even more so for the method of moments since MLE was shown to be far superior in the 
simple power law section of this TP). Thus, an important scientific question is, “For what values of E k will 
these distributions likely begin to overlap for a particular detector?” 

If the concept TSC with its 40-percent Gaussian response function is considered and with an 
observing energy range of 20-5,500 TeV, then equation (20) provides the CRB for each of the three 
spectral parameters as a function of E k . For example, for the case where 6 = (2.8, 3.3, E k ) and calculating 
the CRB for 75 < E k < 400 TeV using equation (20), a 3a curve describing the approximate width of the 
distribution (histogram) of the ML estimate of as a function of E k can be constructed (3.3 minus three 
times the CRB of a 2 ) for values of E k in the 75^00 TeV range. Sketching this curve versus E k and noting 
where it begins to cross the line a, = 2.8 suggests the value of E k where the lower 3a point of the 
distribution of the ML estimate of for this concept TSC will likely begin to overlap that of ct| 
(fig. 26(a)). Also shown in figure 26(a) is the case when the resolution is set to 20-percent and also zero- 
percent (ideal detector), along with three addition dashed curves for the situation where the TSC's 
collecting power is halved. Similar curves are provided in figure 26(b) for the case where (x~, = 3. 1 , giving 
a spectral break size of 0.3. Obviously these figures are no substitute for statistical hypothesis testing and 
furthermore do not consider the two errant estimation possibilities (EE1 and EE2) that also suggest a 
simple power law in favor of a broken power law, but nevertheless still provide important information. 


( a 2 - 3<j crb i Versus E k for TSC — ■ 
and Half-Sized TSC ■ ■ • 



Knee Location - (TeV) 
(a) (0.5 Breaksize) 

Figure 26. 


( a 2 - 3o crb , Versus E k for TSC — 
and Half-Sized TSC ■ ■ « 



Knee Location - E k (TeV) 

(0.3 Breaksize) 

to cross that of a t 
break size is 0.3. 


<b> 

CRB used to estimate where the distribution of a-, begins 
when (a) the spectral break size is 0.5 and (b) the spectral 
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Note, too, that these figures represent a best-case scenario in that they are constructed using the 
CRB calculated from equation (20) which of course is not quite attained in practice, especially for the 
larger values of E k and when the break-size is 0.3 as previously discussed. Thus, the actual overlap point 
would occur sooner (that is, for smaller values of E k ) because the variance of the ML estimator (or any 
other unbiased estimator of ) will be larger than the CRB. 

Furthermore, figure 26 suggests that the 40-percent TSC is roughly equivalent to an ideal detector 
of half its size in terms of measuring a 2 , while similar studies show this 40-percent TSC to be roughly 
equal to a 20-percent resolution detector of half its size in terms of estimating oq and E k .~ Consequently, 
instrument designers should consider first maximizing collecting power and then improving energy 
resolution, whenever possible. 

It is important to realize that while raising E 2 to higher values in this analysis offers no benefit for 
this proposed TSC because of its previously stated collecting power, lowering £j does significantly benefit 
the measurement precision of oq and, because of their correlation but to a lesser degree, E k . However, 
lowering £, offers no benefit in the measurement precision of Ob when using the ideal detector and virtu- 
ally none when using a real detector. These results are presented in table 4 for the case 0 = (2.8, 3.3, 100 
TeV) and in which £j is lowered incrementally from 20 to 1 TeV and once again illustrates the utility of the 
CRB determined by equation (20). Results for the case when the break-size is reduced to 0.3 are similar. 


Table 4. Effect of lowering on the CRB for the TSC-sized detector with 0 , 20-, 
and 40-percent resolution Gaussian response function. The number 
of events N-, above E k is 2,255 for all values of E x . 




CRB-Ideal 

CRB-20% Gaussian 

CRB-40% Gaussian 



Detector (0 

%i 




Detectoi 


(TeV) 

N 


KM 

£*(TeV) 



F*(TeV) 

«i 

a 2 

E„( TeV) 

1 

11,468,838 

H 

Mimjl 

5.98 


0.0578 

7.89 


0.0691 

9.92 

5 

632,364 

EH 


6.09 


0.0579 

8.16 


0.0697 

10.49 

10 

181,152 

EH 


6.25 


0.0581 

8.56 

ill 

0.0704 

11.34 

15 

86,989 

HI 

0.0486 

6.42 


0.0584 

9.03 

0.0126 

0.0713 

12.35 

20 

51,576 

0.0117 

0.0486 

6.60 

0.0144 

0.0586 

9.58 

0.0199 

0.0722 

13.56 


1.4 Analysis of Multiple Independent Data Sets 

The ML theory required to estimate spectral parameters from an arbitrary number of data sets 
produced by science instruments having different observing ranges, different collecting powers, and differ- 
ent energy response functions is developed in this section. Application of this methodology will facilitate 
the interpretation of spectral information from existing data sets produced by astrophysics missions having 
different instrument characteristics and thereby permit the derivation of superior spectral information based 
on the combination of data sets. Furthermore, this procedure is of significant value to future astrophysics 
missions consisting of two or more detectors by allowing instrument developers to optimize each detector's 
design parameters through simulation studies in order to design and build complementary detectors that 
will maximize the precision with which the science objectives may be obtained. 
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This extension of the methodology developed in the previous sections to multiple data sets was 
motivated by such an application and is presented as an example in which two detectors, both assumed to 
have Gaussian response functions but different energy resolutions and observing ranges, were modeled 
separately and then in a collaborative effort to estimate the single parameter of a simple power law energy 
spectrum. A succinct comparison of the benefits from using the sets in concert is measured in terms of 
variance reduction of the estimator, as well as any biases resulting from poor statistics in one or both of the 
individual data sets that may be reduced when considered in combination. 

The ML theory necessary for application to multiple astrophysics data sets is derived here for two 
independent data sets, A and B, produced by instruments having different ( 1 ) observing ranges, (2) collect- 
ing powers, (3) energy response functions, and (4) energy resolutions. These two data sets will be used 
together to estimate the energy spectra information and thereby benefit from the strengths of each detector, 
whereas, singly, they may be inadequate for achieving the science objectives. In practice, the data sets must 
be corrected for systematic errors in the energy response of the instruments in order to achieve the ultimate 
accuracy of the final spectra information based on the combination of astrophysics data sets. Generaliza- 
tion of this approach to more than two independent data sets then follows by induction. 

To extend the ML theory to handle data sets A and B simultaneously, we begin with the probability 

density function for the data set of instrument responses A = { _Vj , a'-i , ■ • • , .v,y } , given bv 

A. * 


Sa<*,-:©) = J g A ( Xj | E, p A ) <p(E\Q)dE , / = !,■-, N A , 
so that the likelihood function is 


( 21 ) 


n a 

Me)=n J g A U;| E,p^)<P(E-Q)dE , 
'=> L*A 


( 22 ) 


where 0 denotes the vector of spectral parameters of an arbitrary energy spectrum, <ME:Q), to be estimated; 
^ A * s number of detected events from observing range, R A , of instrument A having response function, 

g A , and energy resolution p A , so that the corresponding objective function is 


0 A ( 0) = - log [Z. A (0)] (23) 

and the ML estimate 0 A being that value of 0 for which 0^(0) is a minimum, is obtained from equation 
(24) using the simplex search algorithm as 


e A = mmO A (0, . 


(24) 
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The likelihood function and objective function for data set B are similarly defined and because data 
sets A and B are assumed independent, the likelihood function for the two sets considered simultaneously 
is the product 


L ab (Q)= L a (B)L b (Q) , v ' 

so that upon taking the logarithm of each side gives the objective function as 

0AB< e ) = O A( O ) + O B< 0 > (26) 

and the vector of ML estimates based on both data sets considered simultaneously is 

0 AB = minO AB (0) = min[(2 A (0) + C> B (0)] • (27) 

{e} {e} L 

This procedure outlined above for two data sets is readily extended to k independent data sets so 
that any one of the (2 k - 1 ) possible combinations of the data sets acting together provides the ML estimator 
0 5 obtained by applying the Nelder-Mead algorithm to equation (28) as 


e.s 


= min 

{©} 


■ I 0/0) . , 

j € S 


(28) 


where S can be any of the (2 k - 1 ) nonvoid subsets of the set of integers {1,2, ... , k}. 

The statistical properties of the derived ML estimator can then be studied for each simulated 
scenario and, in particular, those relating to (1) bias, (2) variance, and (3) normality can be rigorously 
investigated using graphical procedures and appropriate statistical techniques. 1 *'* — *■ 

The CRB for the estimator of the vector of spectral parameters, 0, using two independent data sets, 
A and B, is derived in appendix B and shown to be the diagonal elements of the inverse of the covariance 
matrix, I, whose l v elements are 


Iij(B)=N a - 


'3 log [g A (.v;0)l x 3 log U A U;0)] 
30/ 30; 


w ,'3 log [,g k (v;0)] . d log [g B (y;0)] 
+ N ‘ { Mi 36 ; 


(29) 


*Some statistical test procedures depend on the outcome of the normality test for estimators. 
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and where ( I ) the notation < . > denotes “expected value” and in practice can be accurately computed using 
Gaussian quadratures; (2) the partial derivatives are also numerically evaluated using equation (15); and 

(3) g A (A-;0) is the probability density function for data set A = {x j J= 1 V A }, given by equation (2 1 ) and 

similarly for data set B = [y jt j = 1, )V B }. As noted in appendix B, the CRB given by equation (29) is 

readily extended to k independent data sets and provides a vital check on the performance of the derived 
ML estimation procedure. If the simulations show the ML estimator of the spectral information to be 
unbiased and also attains the CRB for a given spectrum-instrument combination, then this ML estimator 
will be the best (minimum variance) unbiased estimator possible from combining multiple data sets for that 
particular astrophysics mission scenario. 

Furthermore, when this ML procedure is used in the design phase of an instrument and if the 
simulations show 0 ML is unbiased and attains the CRB for the science mission under consideration, then 
equation (29) can subsequently be used directly to evaluate the relative merits of various instrument design 
parameters (measured in terms of their impact on reducing the statistical error in measuring 0 ) without 
performing additional simulations. An example is given in section 1.6. This is of tremendous practical 
benefit because equation (29) can be evaluated in mere seconds while the equivalent information based on 
Monte Carlo simulations can take several days to obtain. Extending the ML procedure to multiple data 
sets, where in practice each may contain 10^ to lCf* events for a given science mission, is quite challenging. 
Practical studies conducted in Howell 1,2 typically required >1,000 simulated missions to obtain statisti- 
cally meaningful inferences about a single detector design parameter under study. Such a simulation run 
would last >12 h when the broken power law spectrum was assumed (considerably less for the simple 
power law), largely because of the vast number of numerical integrations required in evaluating the 
objective functions and the many steps taken by the Nelder-Mead search for 0 ML . 

1.5 Example of Estimating a, of a Simple Power Law Using Two Data Sets 

The benefits of having the additional nuclei (Z > 8) from a proposed thin sampling calorimeter 
(TSC) for ACCESS at energies above the transition radiation detector (TRD) energy range was recently 
investigated 14 using this methodology and serves to illustrate the application and utility of the approach. 
Specifically, the value of having these calorimeter-provided nuclei in addition to the expected number of 
TRD events was investigated, and the approach was to measure the benefit of including these additional 
events in the estimation of the single parameter, a v of an assumed simple power law as compared to not 
including them. 


In this example the TRD geometry factor is assumed to be 7.58 m 2 -sr with energy response saturat- 
ing at 20 TeV/n and for an observing period of 1 ,000 d. Furthermore, it was decided to consider the species 
Ne-S with an assumed differential spectral index of 2.39, along with a calorimeter geometry factor of 
5 m--sr, providing 593 events >20 TeV/n on average for the calorimeter for a 1,000-d observing period. 

Additionally, a Gaussian response function was assumed for both detectors, and the TRD was 
assumed to have the same linear mean response as the calorimeter but with a constant 35-percent energy 
resolution over its observing range £j to 20 TeV/n, while the TSC was assumed to have a constant resolu- 
tion of 40-percent over its nonoverlapping observing range 20 to 2,000 TeV/n. This upper limit of 2,000 
TeV was chosen for the calorimeter because the number of detected events >2,000 TeV is negligible for a 
detector of this collecting power. Then, E ] is taken to be different values ranging from 0.5 to 7 TeV as part 
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of a parametric study in this illustration. Of course, the upper energy boundary of 20 TeV/n could be 
extended to a higher value if the saturation limit of the TRD could be improved. 

In the simulation, GCR event energies, £ ( -, are simulated from the TRD observing range, £t RB > 
= [£,, 20 TeV], and where the number of these events, jV trd , is a function of the observing range, geom- 
etry factor, and observing time. Then, for each of these simulated incident energies, Ej, a response, Xj, is 
simulated according to the assumed Gaussian response function and with 35-percent energy resolution, 
pj RD . The calorimeter events are simulated in a similar manner using its defining parameters, /V Cal , R Cal , 
and p Cal , and the Gaussian response function. 

Performing the simulation once defines a so-called astrophysics “mission” and provides a single 
ML estimate of oq for each instrument by solving equation (24) with 0 = oq and all instrument characteris- 
tics appropriately modeled through the formulation of the detector response function, g, for each detector, 
and then solving equation (27) when the two data sets are used in combination in the ML estimation 
procedure. Notationally, a Cal and a TRD are the ML estimates of oq when the calorimeter and TRD are 
considered as standalone instruments, respectively, and a Both is the ML estimate of a x when the two 
instruments are used in combination. The simulation is repeated for many missions, each time generating 
a new set of incident energies from the assumed simple power law spectrum for each instrument and their 
respective simulated energy deposits according to the detectors’ response functions, estimating a, using 
each data set separately and then in tandem. The statistical behavior of these estimates relative to plausible 
design variation of response function parameters is then studied. 

Figure 27a shows the ML estimates a TRD , « Cal , and a Both of oq for 25 simulated missions in which 
£, was first set to 5 TeV/n, thus providing 5,275 events on average for the TRD. Figure 27b shows the 
effect of lowering £, to 3 TeV/n, providing an average 1 1 ,660 events for the TRD, and then lowering £, to 
1 TeV/n, providing an average of 56,920 events for the TRD for 25 missions (fig. 27c). We observe that as 
£j is lowered, thereby increasing the number of TRD events as well as extending its effective observing 
range, and hence “lever arm” effect in the estimation of cq, the calorimeter has an ever-diminishing role in 
its contribution to as illustrated by the and 0t Bot j, estimates nearly coinciding in figuie _7c. 

♦ — ^TRD ® a Both ^ ^Cal 



Figure 27a. ML estimates of a, for 25 missions, with £, = 5 TeV (5,275 events for TRD 
and 593 for calorimeter). 


31 



ML Estimate 



Figure 27c. ML estimates of a, for 25 missions, with £, = 1 TeV (56,920 events for TRD 

and 593 for calorimeter). Note that and a Both are virtually indistinguishable 



The simulation was run for 1,000 missions for each of several values of the E x parameter with 
summary statistics provided in table 5 for the calorimeter acting alone; table 6 for the TRD acting alone, 
and table 7 for the two data sets acting in concert to estimate a,. Tables 5 and 6 show that the calorimeter 
and TRD as standalone instruments each provide unbiased spectral information (recall a, = 2.39 for this 
study) and attainment of the CRB. 


Table 5. Summary statistics of a Cal based on 1,000 simulated missions of the calorimeter. 


f, (TeV) 

E 2 (TeV) 

Calorimeter 
Events N ( E > £,) 

Mean 

Standard 
Deviation a Cal 

CRB 

20 

2000 

593 

2.39 

0.066 3 

0.067 a 


a The appearance that o Ca| is less than the CRB is an artifact of the finite number of missions 
in the simulations. 


Table 6. Summary statistics of cty RD based on 1 ,000 simulated missions for each value 
of £j of the TRD alone (all events between E { and 20 TeV). 


E , (TeV) 

Saturation 

(TeV) 

TRD Events 

Mean 

Standard 
Deviation o 1RD 


7 

20 

2968 

2.39 



5.2 

20 

4947 

2.39 

0.068 

0.068 

5 

20 

5274 

2.39 

0.063 

0.063 

3 

20 

11,658 

2.39 

0.031 


1 

20 

56,919 

2.39 

0.011 

0.011 

0.5 

20 

150,630 

2.39 

0.007 

0.007 


Table 7. Summary statistics of a Both based on 1,000 simulated missions for each value 
of E x of the TRD and calorimeter acting in combination. 


CRB for TRD and Calorimeter (Includes 593 Events >20 TeV) 

TRD 

Calorimeter 

Total Combined 
Events 

Mean 

Standard 
Deviation c Bo(h 

CRB 

f, (TeV) 

Saturation 

(TeV) 

f, (TeV) 

Ej (TeV) 

7.0 

20 

20 

2,000 

3,560 

2.39 

0.058 3 

0.059 

5.2 

20 

20 

2,000 

5,539 

2.39 

0.048 

0.048 

5.0 

20 

20 

2,000 

5,866 

2.39 

0.045 3 

0.046 

3.0 

20 

20 

2,000 

12,249 

2.39 

0.027 3 

0.028 

1.0 

20 

20 

2,000 

57,511 

2.39 

0.011 

0.011 

0.5 

20 

20 

2,000 

151,220 

2.39 

0.007 

0.007 


3 The appearance that a Cal is slightly less than the CRB is an artifact of the finite number of missions 
in the simulations. 






























Table 7 shows that for all scenarios considered, the two data sets acting together likewise provide 
an unbiased estimate of the spectral parameter, a,, and attainment of the CRB computed from 
equation (29) with 0= a x to give 


CRB 


Both 


N< 


TRD\ 


d [gjRpl£^l 

da , 




■ + N, 


Cal\ 


d lOg [.?Cal<- V;g l )1 
da, 
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as the CRB for the TRD and calorimeter acting in combination. The CRB in table 5 corresponding to the 
calorimeter acting alone is obtained from equation (14) and likewise for the TRD acting alone. Numerical 
evaluation of equation (30) for the different scenarios of interest presented in table 7 indicates the two 
detectors acting in concert do indeed attain the CRB for all values of £,. Achievement of this bound, 
coupled with the fact that the derived ML procedure provides an unbiased estimate of the spectral param- 
eter, a h implies that this procedure is the best (minimum variance) unbiased estimation technique 
possible. 

Additionally, histograms of the 1,000-mission ML estimates for each instrument considered 
separately and also in tandem show the distribution of the ML estimator to be approximately normally 
distributed. 


Comparing standard deviations of the different scenarios in tables 5-7 clearly shows the impor- 
tance of not only the number of events or so-called “statistics” but also the observing range of the instru- 
ments; i.e., the “lever arm effect.” The very special case where E } = 5.2 TeV, found by trial and error, shows 
the calorimeter acting alone with only 593 events but with a much greater observing range (20 TeV/n to an 
average maximum event energy of 1,100 TeV/n) is effectively as good as the TRD in terms of estimating 
a \ ( a Cal = 0 067, and o TRD = 0.068 when £j = 5.2 TeV). For this case where £j = 5.2 TeV, the TRD had an 
average ol 4,947 events but with a much shorter observing range. Furthermore, the TRD's energy resolu- 
tion was assumed to be 35 percent, whereas the calorimeter was assumed to be worse with a constant 
40-percent energy resolution. We also note the fairly intuitive result in table 7 that a Both = 0.048 which is 
-0.068/ V2 when the calorimeter and TRD are effectively equal to each other and so cr Both scales roughly 
by V2 for this case. 

Nevertheless, as E x is successively lowered, we note the overpowering impact of the TRD's 
increasing number ot events and growing observing range, and in fact, for E x < 1 TeV, the calorimeter's 
contribution to the variance reduction of the ML estimate of a Bot j 1 is virtually nil. Thus, we see that "how 
low” the TRD can observe plays a major role in assessing the value of these additional calorimeter events 
in the estimation of the spectral parameter, a v It was concluded that since the proposed TRD would easily 
observe to <4 TeV/n, the contribution of the calorimeter’s additional events would be insignificant for 
measuring the spectra for Z > 8. Thus, this study resulted in an important design consideration for the 
charge particle identification detector for the calorimeter and underscores the utility of this ML theory for 
multiple independent data sets in the design phase of new science instruments. 



1.6 Using the CRB to Explore Instrument Design Parameters When Estimating 
0= ( CC| , CLj, E k ) of a Broken Power Law Using Three Data Sets 

The benefit of using three independent data sets in concert to measure the spectral parameters of an 
assumed broken power law proton spectrum will be investigated, and, for specificity, assume the spectral 
knee location is E k = 175 TeV and the slope parameters a, lob below/above the knee are 2.8 and 3.2, 
respectively, to give a 0.4 spectral break-size. 

Next, assume the three hypothetical data sets A, B, and C were collected by instruments having 
observing ranges, collecting powers, and stochastic response functions defined in table 8 so that data set A 
provides information regarding only the spectral parameter (X\ \ B about Ofj, o^, and E k , and C, only about 
oc~i. It is further assumed that each data set was produced by a detector having the same lineal mean 
response as the TSC previously introduced in section 1 . 1 so that their response functions may be visually 
compared in figure 28. 


Table 8. Data sets and associated response functions, with CRB for all possible 
combinations of data sets. 


Data Set Combinations and 
Instrument Response Function 

Observing 
Range (TeV) 

Average 
Number 
of Events 

Cramer-Rao Bound lor 
a, a, E t (TeV) 

A (Gaussian) 

[1,20] 

150,000 

0.0095 

- 

- 

B (Gamma) 

[35, 5,000] 

44.000 (2,000)* 

0.0225 

0.0780 

34.24 

C (Broken Gaussian) 

[500, 10,000] 

500 

- 

0.1050 

- 

AB 



0.0087 

0.0744 

26.13 

BC 



0.0220 

0.0626 

30.23 

ABC 



0.0087 

0.0607 

23.05 


*44,000 events above 35 TeV of which 2.000 are above £*=175 TeV. 


In this scenario, data set A was produced by an instrument having a Gaussian response function and 
with constant 40-percent energy resolution; B by an instrument having a gamma response function and 
thus capable of describing a wide variety of shapes with right-hand skewness (outer curve from the right in 
fig. 28); and C was measured by a detector having a “broken Gaussian” response function consisting of 
two blended normal distributions (middle curve from right in fig. 28) suggested by J. Ormes, 2000, private 
communication, for its closeness to the Gaussian response function but with a tail region as desired. These 
latter two response functions were introduced to address the possibility that some detector response func- 
tions may exhibit a right-hand skewness and might add to the difficulty of the ML estimation task. Note 
that while the gamma response function used here also has a constant energy resolution of 40 percent, the 
broken Gaussian actually has a 4 1 -percent resolution because of the added skewness while keeping the rest 
of the distribution matching the Gaussian. 

Furthermore, in the construction of these hypothetical data sets, the number of events in each data 
set was chosen so that each data set alone provides unbiased spectra information and approximate 
attainment of the CRB, 1 * 2 as discussed in sections 1.1 and 1.2, so that any combination of the data sets 
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Figure 28. Gamma, Gaussian, and broken Gaussian detector response functions to 40-TeV proton. 


should also attain the CRB. Consequently, the CRB can be calculated from equation (56) for the 
combination of data sets and used directly to investigate the relative merits of instrument design parameters 
without performing the simulations. Note that while A and C provide no information about the knee loca- 
tion E k , they do provide an improvement in the measurement precision of E k when they work together with 
data set B. 

To illustrate the utility of the CRB for multiple independent data sets regarding design parameter 
studies, suppose the collecting power of instrument C is increased by a factor of 10 with the results given 
in table 9 which shows the merits of collecting power on the measurement precision of the spectral param- 
eters for the combinations of data sets C, BC, and ABC. While C only contains information about ct?, the 
measurement precision (standard deviation) of E k reduces from 30.23 to 23.6 TeV for the BC combination 
and from 23.05 to 17.47 TeV for the ABC combination. 

Table 9. Number of events from table 8 in data set C is increased by a factor of 10. 


Instrument Combination 

Observing 
Range (TeV) 

Average 
Number 
of Events 

Cramer-Rao Bound for 
a, a, E. (TeV) 

A (Gaussian) 

[1,20] 

150,000 

0.0095 

- 

- 

B (Gamma) 

[35, 5,000] 

44,000 (2,000)* 

0.0225 

0.0780 

34.24 

C (Broken Gaussian) 

[500, 10,000] 

5,000 

- 

0.0332 

- 

AB 






1 BC 




0.0306 

23.60 

ABC 







*44.000 events above 35 TeV of which 2,000 are above 5^=175 TeV. 
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Next, if the number of events in data set B is reduced by a factor of 2 and those in C increased by a 
factor of 10, we can once again explore the impact of various collecting power options shown in table 10. 
Finally, if we keep the collecting powers the same (B reduced by a factor of 2, C increased by a factor of 
10), and additionally improve the resolution of the instruments for A and C to 30- and 35-percent 
resolution, respectively, we see the impact of detector energy resolution on measurement precision of the 
spectral parameters, as shown in table 1 1 . 


Table 10. Number of events from table 8 in data set B is reduced by a factor 
of 2 and those in C increased by a factor of 10. 


Instrument Combination 

Observing 
Range (TeV) 

Average 
Number 
of Events 

Cramer-Rao Bound tor 
a, a, E k (TeV) 

A (Gaussian) 

[1,20] 

150,000 

0.0095 


- 

B (Gamma) 

[35, 5,000] 

22,000(1,000)* 

0.0318 

0.1104 

48.42 

C (Broken Gaussian) 

[500, 10,000] 

5,000 

- 

0.0332 

- 

AB 



0.0091 

0.1048 

35.86 

BC 



0.0302 

0.0318 

31.81 

' ABC 



0.0090 

0.0317 

22.36 


‘22,000 events above 35 TeV of which 2,000 are above £ t = 1 75 TeV. 


Table 1 1 . Number of events from table 8 in data set B is reduced by a factor 

of 2 and those in C increased by a factor of 10, resolution of A improved 
to 30 percent and C to 35 percent. 


Instrument Combination 

Observing 
Range (TeV) 

Average 
Number 
of Events 

Cramer-Rao Bound for 
a, a t E k (TeV) 

A (Gaussian, 30%) 

[1,20] 

150,000 

0.0083 

~ 

- 

B (Gamma. 40%) 

[35, 5,000] 

22,000 (1,000)* 

0.0318 

0.1104 

48.42 

C (Broken Gaussian, 35%) 

[500, 10,000] 

5,000 

- 

0.0321 

- 

AB 



0.0080 

0.1047 

35.58 

BC 



0.0301 

0.0309 

31.70 

ABC 



0.0080 

0.0307 

22.02 


*22,000 events above 35 TeV of which 2,000 are above F„=175 TeV. 


Obviously, the number of possible parametric studies are numerous, but the preceding sample 
investigation illustrates the utility of this procedure and the derived CRB for multiple independent data sets 
when considering the design of new, complementary detectors. Furthermore, the CRB in tables 8-1 1 were 
computed using equation (56) in <1 min, whereas the equivalent information based on Monte Carlo simu- 
lations would take > l week to obtain. However, when the CRB is not attained in practice, simulations must 
be used to estimate the real performance benefits when using multiple independent data sets. 
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2. CONCLUSIONS 


This TP investigates the statistical properties of the ML estimator of the single parameter of a 
simple power law energy spectrum and presents the conditions under which this estimator attains the three 
desirable properties: (PI ) consistency, (P2) efficiency, and (P3) asymptotic normality. A comparison with 
the estimation procedure called the method of moments is also included and shows the ML estimator to be 
far superior in terms of these three desirable properties. 

The properties of the ML estimator of the broken power law energy spectral parameters and the 
conditions under which P 1 , P2, and P3 are attained are also investigated under a wide range of parametric 
values. A crucial result of this research and necessity for investigating P2 is the derivation of a closed-form 
expression for the CRB of the broken power law distribution presented in appendix A. Another critical 
result is the calculation of the CRB using equation (20) for the broken power law distribution and equation 
( 14) for the simple power law' distribution when events are measured by a real detector having response 
function, g , and energy resolution, p. While this study considered an instrument having a Gaussian re- 
sponse function and with resolutions ranging from zero to 50 percent, any response function, g, can be used 
in equations (14) and (20) to calculate the CRB, such as various skewed distributions and others having 
nonconstant energy resolution as illustrated in section 1.6 and also discussed in references I and 2. Much 
insight into the estimation task of power law spectra information can be gleaned from the CRB as illus- 
trated in this TP, as well as the fact that it provides a stopping rule in the search for the best (minimum 
variance) unbiased estimator of power law spectra information. 

Simulations were conducted in parallel with these analytical results and the CRB played an unex- 
pected but valuable oversight role by signaling errors in the simulation code when the simulations occa- 
sionally produced the impossible result of the ML estimator having a standard deviation smaller than the 
CRB. Furthermore, it is likely that these simulation errors, while small in practical terms, would have gone 
unnoticed if not for having the CRB available for comparison. 

Additionally, several detector design parameter studies are included in this research and it is hoped 
that those designing instruments to measure power law spectra information will benefit from these studies. 
Additionally, this analysis should benefit those wishing to apply these techniques to the estimation of 
spectra information from existing data sets, which requires a modified likelihood function to handle the 
realistic situation in which the range of integration [£j, E 2 ] in the objective function is unknown, and is 
discussed in detail with examples in references 1 and 2. 

The MLE procedure and companion analytical techniques to estimate spectra information from an 
arbitrary number of astrophysics data sets produced by vastly different science instruments is presented 
and demonstrates how complementary astrophysics missions can work in concert to achieve science goals. 
Additionally, the CRB for an unbiased estimator of spectra information based on these multiple indepen- 
dent data sets is derived in appendix B and provides a means of assessing the accuracy of different estima- 
tion techniques and, furthermore, provides a stopping rule in the search for estimation methodologies 
when the CRB is attained in practice. Several examples illustrating this ML method and utility of the CRB 
for multiple independent data sets are included. 
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APPENDIX A— CLOSED-FORM EXPRESSION FOR THE CRAMER-RAO BOUND 
OF THE BROKEN POWER LAW 


A closed-form expression for the information matrix 1(0) of the broken power law distribution is 
derived, and the inverse I -1 (6) is defined as the CRB. This bound corresponds to the so-called ideal 
cosmic-ray detector having perfect energy resolution and has tremendous utility because it sets the limit on 
the precision with which any conceivable detector of equal collecting power can measure the three broken 
power law spectral paramaters. Furthermore, it provides a means to tune the integration and differentiation 
parameters in the numerical algorithm for evaluating equation (20) for real detectors because the CRB 
determined from equation (42) is exact and must equal equation (20) when the detector’s energy resolution 
p->0. 


A.l Derivation 


The broken power law probability density function for GCR event energy, E , is given by 
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over an energy range [Ej, E~,\ that does not depend on the spectral parameters 0 — (oq, Qb, E ), and the 
normalizing coefficient A(0) is: 
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and writing \og[<p B (E;d)] as 
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where 5, and S 2 are indicator functions defined in equation (37), gives 
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Taking the derivative of equation (35) with respect to a, and a 2 is straightforward and gives 
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where < . > denotes mathematical expectation and the abbreviations 
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are introduced. Taking the derivative of equation (35) with respect to E k and using Leibnitz's rule (see 
Professor Bierens' "check” of equation (39) in section A.2): 
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and the three partial derivatives of log[A(0)] defined in equations (37) and (39) are 
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As a check, the expectation of the partial derivatives of (j)g with respect to the parameters oq, os. 
and E k defined in equations (36) and (39) are indeed seen to be zero when the derivatives A,, A 2 , and A k in 
equation (40) are used in conjunction with the expected value terms defined in equation (41), along with 
(5|) = Pr{£< £ /t },and(5 2 ) = l — Pr{£< E k }, and 
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Finally, the matrix elements defined in equation (42) are multiplied by N to give the information 
matrix 1(0). It is interesting to note that 1(0) contains the intuitive terms N<8 { > and N<8 2 >, which are the 
expected number of events below and above the spectrical knee E k , respectively. 

This completes the derivation of the closed-form expression of the information matrix 1(0) corre- 
sponding to the so-called ideal detector having perfect energy resolution. Its inverse I _l (0) gives the CRB, 
and the square root of the diagonal elements of 1“ 1 (0) is the CRB for the three spectral parameters. In 
practice, one should verify that the difference of the covariance matrix C and the CRB is positive definite 
by checking that the eigenvalues of the matrix [C- I -l (0)] are all positive. If the eigenvalues are zero, the 
CRB has been reached and the ML estimator has attained property P2 (efficiency) by achieving the CRB. 

The computer program MATHEMATICAL was used to provide equations (40H43) and performed 
the important check 


0 = (/\| - L e 8\ j 
< 0 = i^A 2 — 


that confirms the correctness of the terms leading to the construction of the information matrix 1(0). Fur- 
thermore, a format feature of MATHEMATICA called FortranForm was used to write these equations in 
FORTRAN source code for implementation in the overall simulation program and thereby eliminated the 
possibility of introducing human error in transferring the equations into the computer program. 

A.2 Check on Equation (39) by Professor Bierens 16 

“Equation (39) is correct if (38) is correct, which in its turn is the case if (38) holds for fixed p and 
q. In the latter case (38) states that the derivative operator d/da can be moved inside the integral, which 
requires the application of the dominated convergence theorem. The conditions for the latter are: 

1 . For each x in [p,q ], fix, a) is differentiable in a , possibly except for x in a subset with Lebesgue 
measure zero. 

2. There exists a function b(x,a), say, with finite integral over [p,q] such that sup{n[f(x,a+\/n) -fix,a)\ 
< b(x,a) and sup{«|/fx,a-l//i) -fix,a)\ < b(x,a), where the “sup” is taken over all n> 0. 

In (39) you have applied (38) for the cases q = a and p fixed, and q fixed and p = a, so that (39) is 
correct if the conditions 1 and 2 hold for a — q and a = p. In the case q — a and p fixed, the log-density fix,q ) 
is differentiable in q for x < q, but not for x = q. Similarly, in the case q fixed and p = a,fix,p) is differen- 
tiable in p for x > p but not for x = p. Thus condition 1 holds, with {p,q } the subset with Lebesgue measure 
zero. Moreover, since the left and right partial derivatives of fix,a) to "a" are bounded, condition 2 holds as 
well. Consequently, (39) is correct." 
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APPENDIX B — CRAMER-RAO BOUND FOR MULTIPLE INDEPENDENT DATA SETS 


Derivation of the CRB using two independent data sets, A and B, considered simultaneously to 
estimate a single spectral parameter, 0, follow along the lines of the proof of the CR inequality in Hogg and 
Craig 17 for a single data set. Starting with the probability density function for detector A's response (e.g., 
energy deposit) data set A = {jc i -|i = 1,---,^V a } as 

g A (x r ,e)= j g A ( Xi |£,p A )fl£;0)</£\ /= 1, -■•, JV A (45) 

*A 

and similarly for data set B = jvy| j = J as 
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and also for data set B, their product gives the joint probability density function of the two data sets 
considered simultaneously. Assuming the derivative operator dJdd can be moved inside the integral, 
differentiating their product with respect to Ogives 
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Next, define the random variable, Z, as 


(48) 
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It follows from equation (48) that the mean of Z is zero, or using expected value notation, <Z> = 0. 
Moreover, Z is the sum of (N A + N B ) mutually stochastically independent random variables, each with 
mean zero and consequently with variance* 
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Next, assume there exists an unbiased statistic, (7, of the parameter 0 that is a function of the (N A + N B ) 
instrument responses comprising data sets A and B considered acting together, so that 
U = u (-Y| , • • • , x Na VAf B ) and hence 
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and differentiating both sides with respect to 6 gives 
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*The expectation of the cross-product terms is zero because of independence. 
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This shows that <U Z> = 1 , and, because <UZ> = <U><Z> + r Oj/ cr z , where t is the correlation coefficient 
of U and Z, we have 


1 — 0 x 0 + TOfjOz or r= 

G u°z 


(53) 


Since x~< 1, it follows that oq > 1 !a\ and therefore 
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which establishes the CRB for two independent data sets, A and B, for an unbiased estimator of the single 
spectral parameter of an energy spectrum. This bound readily extends to k independent data sets by induc- 
tion and furthermore confirms that “it always pays to use additional data sets” because of the additional 
variance reduction and hence improvements in measurement precision, as long as the additional data sets 
provide unbiased information. 

The above derivation generalizes to the case where 0 is a vector of spectral parameters by the 
additivity of information matrices 1 ^ so that the CRB for the individual parameters comprising 0are the 
diagonal elements of the inverse of the information matrix, I, having ly-elements 
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and, as in the single parameter case, equation (55) is readily extended to k independent data sets. 
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For the interesting case where a detector does not observe events in an energy range represented by 
one or more of the spectral parameters, the partial derivatives with respect to those parameters are zero. For 
example, suppose detector A observes only events below the spectral knee, detector C only observes events 
above the knee, and detector B observes both above and below the knee. Then the matrix to be inverted to 
obtain the CRB for an unbiased estimator of 0 = (a,, o^, E k ) is 


I ij= N A 
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where the notation in the second term of the right-hand side of equation (56) defines 0] = Ofj , 0 2 = a 2' anc * 
=E k as before and where the integration range for this term must be split at E k in its numerical evaluation. 
An example usin^ equation (56) is provided in section 1.6 of this TP. 
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